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Abstract 

In this paper we present algorithms for approximating real band- 
limited signals by multiple Gaussian Chirps. These algorithms do not 
rely on matching pursuit ideas. They are hierarchial and, at each 
stage, the number of terms in a given approximation depends only on 
the number of positive- valued maxima and negative- valued minima of 
a signed amplitude function characterizing part of the signal. Like 
the algorithms used in [6j and unlike previous methods, our chirplet 
approximations require neither a complete dictionary of chirps nor 
complicated multi-dimensional searches to obtain suitable choices of 
chirp parameters. 

1 Introduction 

A real-valued signal / : M — M is said to be band-limited if its Fourier 
transform, usually denoted by /, has compact support. This constitutes a 
broad class of signals sometimes referred to as the Paley-Wiener space being 
of paramount importance in the applications as most of human and natural 
phenomena should exclude infinite frequencies. The classical Paley-Wiener 
theorem states that any band-limited signal with finite energy [i.e. such 
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that Jjg dt is finite) is indeed the restriction of an entire function of 

exponential type defined in the whole complex plane [T^. Consequently, it 
would suffice to know a small part of such a signal in order to be able to 
extend it arbitrarily to any open set of the complex plane. However, this is 
unrealistic at the time being (even if some progress has been made recently, 
see for instance [IB] and references therein) and in practice, one has still to 
focus on their analysis. 

As the Fourier transform is supposed to be of compact support, it may 
seem a good idea to express numerically band-limited signals relying of opti- 
mized algorithms like the well-known Fast Fourier Transform (FFT). Espe- 
cially, Shannon's sampling theorem ensures that all the signal's information 
can be recovered from the knowledge of a countable collection of samples. 
However, we shall explain in §2.2 that this idea doesn't lead to economical 
(so-called sparse) representations (consult especially [8| in this direction). 
An alternative can come from the decomposition into chirps which read like 
A{t) exp{j(l){t))] the term chirplets has been coined by Steve Mann [13] in 
an attempt to derive an object more sophisticated than wavelets. Briefly, 
if Fourier transform gives a full resolution of frequencies but a null time- 
resolution, wavelets offer a good compromise according to the uncertainty 
principle but still decompose signals onto horizontal rectangles in the time- 
frequency plane [5]. As chirps are endowed with a time- varying phase func- 
tion t 0(t), decomposing a signal into a superposition of chirps means that 
time-frequency curves are now involved. The most usual chirps are the gaus- 
sian/quadratic ones, for which both ln(y4) and are polynomials of degree 
2: we shall work in this framework hereafter. 

Clearly, since chirps are more complex objects, they need more parameters 
to be defined correctly. Indeed, one can roughly count that 6 parameters are 
needed for the definition of both functions, plus 2 other parameters dealing 
with the time and the mean frequency around which the chirp is localized. 
All in all, we found at the end of §2.3 that a signal being written as the 
sum of p + 1 chirps requires at most 6p + 3 parameters to evaluate. In many 
cases, this constitutes a great improvement with respect to classical Fourier 
techniques. 

The main drawback in the setup of chirp decomposition techniques is the 
weight of the computational effort required: according to the literature, most 
of the algorithms rely on matching pursuit techniques [12] where one first con- 
siders a rather complete dictionary of chirps in order to find iteratively the 
ones matching the signal as best as possible (see [H IH [3 HDl El E] and 
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the recent paper [2] motivated by gravitational waves detection). However, 
in a very recent paper, Greenberg and co-authors [S] proposed a completely 
different approach where one gets rid of the dictionary and constructs the 
chirp approximation by means of a simple and easy-to-implement procedure. 
Loosely speaking, given a complex signal in polar form uj i— > a{uj) exp(j(y9(u;)), 
it amounts to seek local maxima of a (we call them for instance uj^) and 
approximate both ln(a) and if by polynomials of degree 2 admitting an ex- 
tremum around uOn- The procedure is first applied to tui, the first maximum 
point of a, in order to derive 2 polynomials Ai and 0i, and then it applies 
to the remaining signal a(a;) exp(j(^(a;)) — y4i(ci;) exp(j0i(a;)) until residues 
become low enough. In the present paper, we propose a more elaborate algo- 
rithm than the one contained inside [B]; it will be presented in detail within 
§3. First, §3.1 will deal with a pointwise selection procedure to compute 
chirps' parameters. Then in §3.2, a mean squares procedure will be pro- 
posed. Both approaches will be tested on an academic example in §4 which 
presents a numerical validation of these algorithms. We insist on the fact that 
the approach is numerically efficient and computationally cheap. Finally, §5 
will be devoted to more "real-life" experiments: trying to decompose a signal 
slightly corrupted by white noise and seeking a chirp inside a stock market 
index. 

2 Preliminaries 

2.1 Specificities of real, band-limited signals 

Our interest lies in an efficient, approximate representation of real- valued, 
band-limited signals t f(t). If /(■) is such a time-dependent signal, it may 
be written in a very general way as 

1 

Vt e R, fit) = - / e^'"^H{uj)du, (1) 

where j = is the unit of the imaginary axis, u stands from now on for 

the Fourier dual variable and the function Ti, rewrites like: 

Vcu G i-n, n), n{uj) = h^{uj) - jHo{uj). 

The corresponding even and odd components of 7i, denoted by Hf.{-) and 
Ho{-), are smooth real- valued functions satisfying 
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^ue{-Q,Q), He{-uj) = He{uj) and Ho{-uj) = -Ho{uj). 

By definition of band-limited, botfi and Ho vanisli identically outside of 
the interval (— moreover, they are assumed to satisfy 

\mv HeiVL - e) = \im HoiVL - e) = (2) 
The symmetries of He{-) and Ho{-) imply that 

Vcj G (-^], fi), = 7^(-cj), (3) 

which is a well-known property of the Fourier transform for real signals (the 
overbar stands hereafter for complex conjugate). The function can be 
recovered from the signal via the classical Fourier transform: 



n{uj) = / e-^'^'f{t)dt. (4) 

J ~oo 

This last identify further implies that for all u, 

He{uj) = 2 cos{ujt)fe{t)dt; Ho{uj) = 2 sm{ujt)fo{t)dt (5) 
Jo Jo 

hold for /e and fo, the even/odd parts of the signal / defined for all t G M 
as follows: 



fe{t) =^l{m + f{-t)) = U-t) and m \[f{t)-f{-t)) = -U-t). 

We regard (jll)-(l5]) as a sanity check on how well we are doing in synthesizing 
/(•) from 7i(-); that is if we employ some algorithm to compute ([1]), approx- 
imately, we then use that computed /(■) to approximately evaluate ([3]) and 
see how well the result agrees with the original function 7i(-). 
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2.2 Standard representations of band-limited signals 



Standard representations of band-limited signals may be obtained by dis- 
cretizing the integral defined in ([1]). If we introduce a discrete (and finite) 
set of frequencies Uk 

kVl 

o^fc = — ,ke{-N,-N + l,...,N}, 

and exploit (Ej), we find that the trapezoidal rule, applied to ([H), yields the 
approximate band-limited function, fN{-), whose value reads 

k=~N+l ^ ^ 

The function is real-valued and satisfies — /Ar(t)| = 0(1/A^^) provided 
the functions -ffe(') and -f^o(') are on [—^2,^2]. More interestingly, /Ar(-) is 
periodic with period, T/v = and is completely determined by its values 
at times = ~^ < n < N — 1 (this is the classical Shannon sampling 
theorem). To see this we note that 




which yields 



n=-N k=~N+l ^ ' \n=-N 



gJ AT 



Now, one observes the following property of the exponentials: 
N-i r , if k^p 

E. {k-p)mT I 
gj JV = / 

n=-N [ 2N , if k = p. 

The identities (Q-iD imply that for -A^ + 1 < p < - 1 

N-l 



^ ^ n=-N 
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while for p = ±N, they yield: 



0- (=f)^Ee-/.(^). (7) 

n=-N n=-N 

This identity ([7j) shows that if we extend to the indices p = ±N, 
then the extension is consistent with the constraint ([2]). This last set of 
identities give an alternative means of computing 7i(-) at the lattice points 
— < P < iV- For completeness, we record relevant identities for the 
coefficients {H^ {^) , Ho (^)), < p < - 1. They are: 



(f ) - - S t (^) (i) , . P S A' - 1, (8) 

^ / n=l 
^ / n=l 

where 



'nn\ def 1 



/iv( 



N,e 



-nn 



, l<n<N, 



and 



(i) K^" (i) - (^)) (^) • ^ ^ " ^ 

The periodicity of /Ar(-) guarantees that /7v,e ('^) = /a^ (~^) 
that /at o ('^) = 0. Moreover, if we extend to p = N, then ([7]) guarantees 
that the extension is consistent with Similarly, if we extend to 

p = N, we find the extension is also consistent with ([2]). The standard 
approach outlined above will, in the limit as N ^ oo, yield the desired 
signal /(■) but is computationally intensive and requires a significant amount 
of data. Our goal here is a more economical approximate representation of 
/(•) which, in many circumstances, requires substantially less data. 
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2.3 Chirplet representation of band-limited signals 

We first note that if He and Ho are defined as in §2.1 and satisfy ([2]), tfien 
7i(-) can be written in polar form (at least at points where it doesn't vanish): 

where the "amplitude" is real, nonnegative and satisfies: 

Aiu;)'=' {H!ico) + H!iu;)Y^' = Ai-co), 
Furthermore, A{uj) is identically zero on |ti;| > Q and satisfies lim A{Q — e) = 
The "phase" u ^ is a smooth, odd function satisfying 

cosip[uj) = and smip[uj) 



A{uj) ' A{uj) 

As our canonical model for A{-) we assume that it has exactly p local 
maxima at the following distinct points, 

< Qi < Q2 < ■ ■ ■ < < ^, (10) 

and that each of these maxima is non-degenerate; i.e. A^'^^Qk) < 0, 1 < k < 
p. The fact that A{—uj) = A{uj) guarantees that A{-) also has non-degenerate 
local maxima located at the symmetric set of points: 

-n< -fip < . . . < -r^a < -^^1 < 0. (11) 

The origin, u = 0, will be either a local maximum or a local minimum of A{-). 
In the former case we shall assume that A^'^\0) < 0. Given the structural 
properties of the amplitude A{-) we attempt to approximate it by means of 
functions Ap{-) of the following gaussian form: 

Ap{uj) = aoe ^-o + ^ I e ^"k +e ^"k 1 (12) 
fc=i ^ ^ 

where > 0, o"a; > 0, and 

< CiJi < 0^2 < • • • < t^p- 

When uj = corresponds to a minimum of A{-) we choose ao = 0. Approxi- 
mate ^(■)'s of the form f|T2l) are not band limited but they have tails which 
decay rapidly as ^ 00 and this is adequate for our purposes. 
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In §2] we present two algorithms for choosing the parameters > 0, cr^ > 
0, and the numbers u^s and show how these selection procedures perform 
on various examples. For the time being, we assume the relevant parameters 
are known and instead of working with 7i(a;) = A{uj)e~^'^^'^\ we replace it 
with 

=^A,(^)e-^'^(-). (13) 

Using this approximation to 7Y(-) we arrive at the following approximation 
for the signal defined by ([1]): 

/•OO 2 



OO 



p 



k=l 



OO 2 



+ e^-^* / e-^+^'^"*-'^^"'=+"»dn (14) 



OO 



/OO 2 
g-5^+j("*-V'(-"'fe+«))^^ 
-OO 

To carry this process further, we replace the terms ilj{±Wk + u) in ( |T4l) by 
their local quadratic Taylor approximations around ±Wk. The result is 

2 2 

tp{uJk + m) = 7fe + tkU H —, p{-uJk + u) = -7fc + tkU —, 

hence 

tk = ipiujk) = ip'i-Uk), (15) 

Kfc = ^'"(cUfc) = -Ip^i-UJk)- 

We note that this last step is equivalent to replacing Hp^i{uj) with 



A;=l 



(16) 



fc=i 



and computing (JT]) with 7i(-) replaced by 7ip,2(-)- The function 7ip,2(') sat- 
isfies 7ip,2(i^) = 'Hp^2{,—^) as required by ([3]) of 7i(-)- The reader will note 
that 7ip,2(") is merely the sum of 2p + 1 complex Gaussian (or quadratic) 
chirps. Exploiting the quadratic approximations defined in f|T5l) when eval- 
uating yields, after a tedious calculation, the approximation /p2(-) to 
/(■): 



/p,2(t) = (27rao)^/^aoe-^ 



1/2^ 2{l + «;Ja^) 



fc=l ^ K k' 



where 



'l+3^^k(Tk) = {l + Kiaiyl'e 



2 2\l/2^2i<^^ 



or equivalently 



1 l^k'^k 

The aforementioned approximation /p,2(') is simply the sum of 1 real 
Gaussian chirps which are characterized by 6p -|- 3 parameters 

• (ao,cro,io), 

• (afc,u;fc,crfc,7fc,4,Kfc), for k = {l,...,p}. 
Finally, 0^ is computed using (|T7j) . 



2.4 The objective of the paper 

Approximations of band-limited signals by real-valued Gaussian chirps of 
the type described in fl2.3p offers efficient representation of radar, seismic, 
and some biological signals. A sampling of the literature on this subject 
may be found in [H [21 13 El El [HI HH [15] and the references contained 
therein. The central problem in establishing such a formula is an efficient 
and natural method for obtaining the parameters characterizing the chirps. 
The "matching-pursuit" and "ridge-pursuit" algorithms have been used by a 
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number of authors; for details see e.g. [H IH El [T21 [IS] . These typically require 
complicated mult i- dimensional searches to obtain the chirp parameters and 
they also demand complete a-priori dictionary of chirps. The procedures 
we advance in ^ to capture the parameters (at, 0"^, tUfc)^^^ require no such 
complete dictionary and rely only on information about the signal amplitude 
A{-). The remaining chirp parameters (7^, tfc? i^k)'k=i obtained from local 
information about the phase, near the points uj^- The algorithms are 
also hierarchical and give our approximations a multi-resolution character. 
Less refined versions of these procedures were advanced in 

3 Parameter Selection for Ap[-) 

In this section we present two different algorithms for the selection of the 
parameters defining Ap(-)^ recall ffT^ . For definiteness, we assume a; = is 
a local minimum of A{-) and thus choose = 0. 

3.1 Pointwise Selection Procedure 

We attempt to choose the parameters so that at the positive local maxima, 
Qi < Q2 < ■ ■ ■ < ^p, one gets: 

Ap{Qj) = A{Qj),A^j^\Qj) = A'-^\Qj) = 0, and A'^\Qj) = A^^\Qj) < 0, 

(3.1), 

for 1 < j < p. Hereafter, we denote by g the Gaussian function centered in 
±uJk: 

g{u;Uk,(Jk) = ie +e J . (ig) 

Then, its first derivative reads 



g'^^\uj,ujk,ak) = - 
and its second derivative is 



-e + -e 

0"fc 



g^'^\uj]UJk,(Tk) 



1 

+ 



{UJ - UJkf 



at 



+ I 



e 
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Hence, for 1 < j <p, the set of equations (3.1)^ rewrites as 



ajg{uj;ujj,aj) ^ A{Qj) -^akg{^j;ujk,crk), (19) 



fc=i 



which boils down to: 



_(tai2! ^ (".7 +"7)' 

ttj-e ^A{Qj)-2_^akg{flj]LJk,cTk)-aje '"'^ . 

k = l 

We differentiate this equahty once, 



•' k=i J 

k^j 

and twice: 



k=l 
k^j 



+aj I ^ - ^^1^ ) e 



0-7 



We have in mind to propose an iterative algorithm to compute a solution to 
the above system of 3 equations in order to get the value of the 3 parame- 
ters; the index n will refer to the successive approximations built out of this 
iterative process. 

• For the first index n = 0, we let 

a? = A{nj), = Q,, and = -^^y for l<j<p. (20) 

• We now suppose the coefficients (cr", 0;", '^j'Yj^i are known; we want to 
compute these parameters at step n + 1. For the first index j = 1, we 
let 
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m+i 



fc=2 



fc=2 ^1 ^ 1 



If > and < 0, then we can compute an auxiliary quantity: 



and cr""^^, cj"'^"'^, and a"^^ are deduced as follows: 



„n+l 



n+1 



n+l 



n+1 



n+l ' 



(21) 



> 0. 



We now assume we have determined (a 



n+l ,n+l _n+l 



0""'''^) for indices 



1 < J ^ Jo < P — 1- We then fix the index jo and computqj at step 
Jo + 1: 



JO 



-«"o+ie 



fc=i 



("io+i+""o+i) 



5^ «fc5(^^jo+i;^fe,^^fc) 

k=jo+2 



then, 



^li jo = P ~ 1, the last sums in ((22|) - ((24|) are not present. 



(22) 
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fe=l 



2 



fc=io+2 

and: 



io 

fc=i 

+".0+1 ^ — 7^^-^^ — e 



p 



If A^^_^i > and C^^i < 0, then the auxihary quantity reads 

7-)n+l »n+l / /in+l 

and compute <^j^„+i and o.'^^i as follows: 



An+1 



j'o + lV jo + l/ io + 1 



, ,n+l 



n+l 
io+1 



^;+\e ^ > 0. 



We thus iterate until satisfactory convergence is obtained. 
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The aforementioned algorithm reahzes a compromise between simphcity and 
efficiency because it stems on inverting only half of the nonlinearities appear- 
ing in the equation f|T9|) and its derivatives in order to keep computations 
easy as it suffices to define the auxihary quantity V to derive fl2Tl) and fl25|) . 
This way of proceeding is related to standard iterative methods in numerical 
linear algebra, like for instance Jacobi or Gauss-Seidel. 

Remark 1 We note that with this method of parameter selection, the fol- 
lowing estimate holds: \A{ijj) — Ap{uj)\ = 0{\uj — in a neighborhood of 
each of the local maxima, Qj, of A{-). This pointwise selection procedure 
generalizes the procedure used in 16^ and captures the behavior of A{-) near 
all maxima simultaneously. 



3.2 Parameter Selection 

Again, we seek to approximate the amplitude A{-) by a function of the form 
ylp(-), recall (|T2|) . and choose the parameters (a^, u;^, ct/c) so as to minimize 
the mean-squares error: 

/oo 
-oo 

/n poo 
A'^{uj)duj-2 A{uj)Ap{uj)duj + / Al{u)duj. 

For definiteness, we again assume that = is a local minimum of A{-) and 
choose ao = 0. We also adopt the short-hand notation 



Gi{u)=\^e +e j = g{uj;LUi,ai), I < i < p, (26) 

and note, for future reference, that for all 1 < i and j < p: 



{Gi, Gj) — {Gj, Gi) 



def 



Gi{uj)Gj{ijj)duj 



2o"jO",-7r 



1/2 



7 



+ e 



' 2{<Ti + <T.) 



(27) 
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A quick calculation shows that \\A{-) — Ap{-)\\'^ is given by 



\\Ai-)-A, 



pU p p 

/ A\co)duj -2j2fkak+Yl ^ 0' (28) 

'^"'^ k=l i,i=l 



where the value of (Gj, Gj) has been given above and 

fk =^ {A, Gk) = r A{Lu)Gk{uj)du;, I < k < p. (29) 



J-Q 

For a given choice of parameters {ui, cri)^=i the ck's which minimize ( 1281) 
satisfy 

Ga = /, (30) 

where Q is the symmetric, positive-definite, p x p matrix whose {i,jy^ entry 
reads: 

Gi,j = {Gi,Gj) = Gj^i, (31) 

ex. is the p X 1 column vector whose i*^ entry is a-j and / is the p x 1 column 
vector whose i*^ entry is /j. The solution to fl5U]) is given by ck = G^^f, and 
this implies: 



A,(-)|P= r A\oo)duj-fG-'f>^. 



Thus, the task before us is to choose (cui, (Tj)f^^ to maximize f^G'^f where 
again G"^ is the symmetric, positive-definite inverse of G defined in (1311) 
above. At first blush, the problem of maximizing Q = f^G^^f looks rather 
daunting but, in fact, this problem has a rather simple structure. 

• For any integer k between 1 and p we let be one of the two parameters 
uJk or (Tfc and we observe that: 



fG-\p,f + 2fG''f,p,. (32) 



Noting that f^G ^ = ol^ and that is the p x 1 column vector 
whose k^^ component is fk^/s,. = {A,Gk,/3^.) and other components zero. 
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we see that the second term on the right-hand side of ([22D is 2akfk,f3k- 
Moreover, if we exploit the identity, 



G = -G ^G,i3k G \ (33) 
we find that the first term on the right-hand side of (132|) is 



2ak^{Gk,Gj) aj-al{Gk,Gk) 



and again the a's are the solution of (|5U]) . These observations imply 
that 



dQ_ 



( 



2a, 



\ 



d{ Gk,G,) 



\ 



Oil {Gk, Gk) 



J 



where again Get = f. Given the particularly simple structure of 
we solve the problem of maximizing Q by a "steepest ascent" type 
algorithm as we explain now. 



So we assume now (a;", cr'^Y^^i being known. With these data, and the 
explicit formula fl26l) for Giiuj), we explicitly compute the (G,, Gj)'s and 
{Gi, Gj)'s at the n^^ data set. Once again /3j = Ui or cTj. We denote 

the results by {Gi,Gj)"' and ^ {Gi,Gj)"' respectively. We sample the 
amplitude A{-) at points 



UJr, 



-N +l<p< N -1, 



and do the same for the functions Gi(-) and ^(O- These are of course 
also evaluated with the n^^ data set and the results are superscripted 
with the index n. We then derive approximate values for fi and |^ 
using the discrete inner products defined below to replace the integrals 
fsee (EHD) 



N+l 



p=~N+l 



pQ 



pQ 
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and 



Next we obtain a" by solving Q'^ol^ = /" in order to compute 



2a? 



dfl 



\ 



3 = 1 



ia:y{G,,G,r , (34) 



and use fl34|) to update the parameters as follows: 



n+l 



duji ' 



Unless 



0, 1 < i < p, Q will increase for A small enough. 



We iterate the aforementioned process till convergence. 



4 A first "academic" numerical example 
4.1 First test 

To see how these selection processes perform we consider the following ex- 
ample. We let 



0, 



-oo < CO < —2 



A{u;) 



(35) 



0, 



2 < CO < oo. 



The origin = is a local minimum of A and the points u = ±1 are the 
maxima of A{-) with v4(±l) = 13.5. The second derivative at these points is 
A^'^\±l) = —36. Our approximations are of the form: 



Ai(uj) = ai { e _|_ g 2<Ti 



(36) 
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The pointwise selection procedure yielded the following results: ai = 13.4515, 
uji = 1.0074 and (Xi = 0.3595, which in turn yields maxAi{uj) = Ai{±l) = 

13.5. The L2 procedure furnished ai = 13.8189, cui = .8974 and ai = 0.2942; 
in this case maxy4i(u;) = Ai{±.9) = 13.8782. The maximal value of Q is 

390.9413 and the norm of A{-) is 395.0055. Figure [H below, shows a 
graph of A{-) along with graphs of Ai{-) for each selection procedure. 

FIGURE 1 : A vs w - black; A1 vs w, Pointwise Procedure - red; A1 vs w, L2 Procedure - red 
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Figure 1: Comparison between pointwise and selection procedures with 
example ( I35l) -( l36|) : original signal A is in black, pointwise algorithm in blue 
and in red. 

The results for this simple example are not atypical of the general case; 
namely, one application of either selection process typically yields an approx- 
imation, Ap{-), which is qualitatively similar to the given amplitude, A{-), 
but may differ quantitatively from A{-). This defect may be overcome by 
repeated applications of the procedure. 
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4.2 Hierarchical Refinements of Selection Procedures 

We assume we have already applied either the pointwise or selection pro- 
cedures n — 1 times and constructed the functions A^^^ki,-), 1 < k < n — 1, 
each which are sums of Nk Gaussians. We let 

ra-l 
fc=l 

where u Ao{uj) is the original amplitude, previously denoted by A{-). As 
defined An{-) is even and rapidly decreasing as \uj\ oo. The principal 
qualitative difference between and Ao{-) is that takes on both 

positive and negative values. For definiteness, we assume now that = is 
a maximum of An{-) and that the typical structure of An{-) is as follows: 

• there are exactly g„ negative- valued, local minima of An{-) at points 

Ka'l < ... <^]^^ 

and each minima is non-degenerate; 

• there are exactly p„ positive- valued, local maxima of An{-) at points 



< f2i < Qr, < . . . < Q„ 
and each maxima is non-degenerate; 

• the minima and maxima are not necessarily interlaced since may 
have positive valued local minima and negative valued local maxima. 

Graphs of Ai{-) = Ao{-) — Ai,o(-) for the two selection procedures used in 
our previous example are shown in Figure [2J Our induction step is to replace 
An{-) by a sum of = g„ + p„ + 1 even Gaussians: 




(37) 
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FIGURE 2: AO - A1 ,0 vs w, Pointwise Procedure - blue, AO - A1 ,0 vs w, L2 Procedure - red 
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Figure 2: Difference Ao{uj)—Aifi{uj): pointwise selection in blue, selection 
in red. 

where the a's, /5's, a's, and a's are positive and the a;'s and aJ's satisfy 

< a;^ < o;^ < . . . < <, < oJ^ < oJ^ < . . . < aJ^„. 

The pointwise selection procedure generates the unknown parameters 
by insisting that AN„,n{-), and match An{-), A'n\-), and 

A^n\-) at the local maxima {0, {^jll^} and local minima {9,1} tU- The 
selection procedure chooses the coefficients to minimize ||74n(-) — 74Ar„^n(-)|p. 
This problem has the same structure as the optimization problem discussed 
in detail earlier, and may be solved by the "steepest-ascent" algorithm. The 
optimal solution satisfies 

< 1 1 A„+i(-) 1 12 = p„(.) - An^A-) 1 1' = I I An(-) I r - Qn,max, (38) 

where (t^z", o:'^; cJ", a") ^ Qn = H^Afn.nlOlP when the a's and /3's are the 
least squares parameters corresponding to the given choice cu", (t", cu", 
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and Qnmax = Hiax Qn- The equality (l38ll implies that 

{uJ.'^ ,g_" ,LJ^ 

n n 

i=o i=o 

and provides us with a stopping criteria for the number of applications of 
the LP' selection procedure. Specifically, we stop when ||yl„+i(-)|p < estop, 
a preassigned tolerance. The stopping criteria for the pointwise selection 
procedure is equally easy; we stop when 

max(v4„+i(u;)) < estop and min(A„+i(cj)) > -estop- (39) 

Ll) Ll) 

Figures [3] and m shows the results of applying the hierarchical L2 procedure 
2,3, and 4 times, respectively. The "black" curves on each figure are the 
original amplitude, A{-), and the "red" curves are the composite hierarchial 
gaussian approximations. In applying the LP procedure twice we went from 
1 to 5 even gaussians; in the third application we went from 5 to 16 even 
gaussians, and in the fourth application from 16 to 47 even gaussians. After 
four applications of the algorithm, the curves become indistinguishable. 



5 Two "real-life" applications 

5.1 Chirp decomposition of a noisy signal 

A band-limited signal slightly corrupted by white noise isn't band limited 
anymore; however, if we set up the preceding algorithms with a value of VL 
being large enough, it may be possible to recover a correct approximation. 
First of all, we set up the following data: consider the even amplitude function 
displaying only two bumps, 

, exp(— alcijP) — exp(— , , , 

A{u) = ^-^ ^; a = 0.8, b = 1.3. 40 

b — a 

We couple it with different phase functions of increasing complexity: first, 
we chosed a cubic one, (/>(a;) = ^ and then (l){uj) = 7r(l — exp(— cj^)) sin(2co'). 
The discretization grid is 512 points from t = —5.1 to t = 5.12. We don't 
insist on the fact that in case the original phase function is a polynomial 
of degree 2, its recovery is exact. 
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A vs w-black and AA vs w, L2 Procedure-red, number of applications of tlie aigoritlim and total number of chirps =2 5 
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Figure 3: hierarchical refinement procedure (to be continued in Fig. H]): 
Original signal A is in black and in red. 

The specificity of these runs is that now, we test the algorithm "blind", 
which means that we just furnish the collection of sample data for both 
amplitude and phase, and we look after the maximum numerically. Similarly, 
the derivatives involved in the parameters calculation are computed with 
finite differences. Corresponding results are displayed in the Figures [5H71 
On the left column, one sees the original amplitudes and phases (in black) 
and their approximations (in blue); on the right one, there are the signal 
as a function of t (in black), its chirplet approximation (in blue) and finally 
the absolute error in log-scale. The recovery of the amplitude ( HOj) by means 
of two Gaussians looks very satisfying and errors are noticeable only in the 
last example where the original signal is corrupted by white noise; the cubic 
phase (Fig. [5]) is approximated in a correct way by second degree polynomials 
around the maxima of A which are well identified. Of course, in case A would 
display several local extrema, a more involved process would have to be set up 
in order to localize them properly inside the vector of samples. Concerning 
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the sinusoidal phase model (Fig. E]), its recovery by means of polynomials 
of degree two is of course rather poor; however, what really matters is that 
it is correct in the vicinity of maximum points of A, and this is just what 
happens (see [IS]). The behavior of these approximations in the t variable is 
good, and the absolute error remains reasonable. 

Fig. [7] deals with a more difficult test-case, namely we set up the ampli- 
tude (HUl) with the sinusoidal phase, we perform the inverse Fourier transform 
to get it as a function of t, we corrupt the resulting signal with white noise 
(which makes it not band-limited anymore). This is the data we furnish to 
the algorithm of §3.1. Obviously, the white noise has to remain small in order 
not to perturb the maxima of A. What we found is that even this require- 
ment isn't enough as the phase function should not be corrupted too much 
in order to let the algorithm approximate it locally by a parabola correctly. 
In this case, the absolute recovery error is bigger compared to the noise-free 
case (see again Fig. [6]). 



5.2 Finding chirp patterns in stock market indices 

The Hang- Seng Composite Enterprise Index is one of the main index for the 
Chinese stock market and the CAC 40 is the leading index for the Parisian 
place. We used the corresponding daily price fluctuations on an interval of 
1024 days (ending on august 29th 2008) in order to check on a real-life case 
whether or not the level of tolerance to noise was acceptable. Clearly, in order 
to avoid as much as possible spurious local extrema in the Fourier spectrum 
of the data, we detrended the prices by a global least-squares interpolatioiJl; 
polynomials of degree 5 and 6 have been used. Moreover, for the CAC 40 only, 
the detrended fluctuation was displaying quite a sharp peak corresponding to 
a periodic cycle. We also used a Basis Pursuit algorithm to remove as much 
as possible the white noise components of these fluctuations, see [3] . On Fig. 
m we display the detrended data (in black), the long-lasting periodic cycle 
of the CAC 40 (in blue) and the chirps we got out of the algorithm of §3.1. 
The case of the Chinese index seem to be quite interesting as the recovery 
seems to be quite sharp. Results on both CAC 40 and HSCEI suggest that 
a qualitative change in the behavior of blue chips quotations occurred after 
the ignition of the so-called subprime crisis/credit crunch (corresponding to 
abscissa t ~ 700). 

^This construction ensures the fluctuation has some amount of vanishing moments. 
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A vs w-black and AA vs w, L2 Procedure-red, number of applications of the aigorithm and totai number of chirps =3 1 6 
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;ure 4: hierarchical refinement procedure, continued from Fig. 
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Figure 7: Pointwise procedure selection on example (HOl) slightly corrupted 
by white noise with sinusoidal phase .^^ 




Figure 8: Chirp patterns in the daily price fluctuations of CAC 40 (top) and 
HSCEI (bottom). 30 



